Comprehensive coverage of prior art in the field of image restoration relevant to the system and method may be found in H. C. Andrews and B. R. Hunt, Digital Image Restoration, Prentice-Hall Signal Processing Series, Prentice-Hall, Inc., Englewood Cliffs, N.J. (1977); W. K. Pratt, Digital Image Processing, 2nd Edition, John Wiley and Sons, NY (1988); H. Stark, Image Recovery, Theory and Application, Academic Press, Inc., Harcourt Brace Jovanovich, Publishers, New York (1987); R. C. Gonzalez and P. Wintz, Digital Image Processing, 2nd Edition, Addison-Wesley Publishing Company, Inc., Advanced Book Program, Reading, Mass. (1987); and R. L. Lagendijk and J. Biemond, Iterative Identification and Restoration of Images, Kluwer International Series in Engineering and Computer Science, Kluwer Academic Publishers, Boston, Mass. (1991). For space-invariant blurs, digital image restoration is associated with the solution of two dimensional ill-posed convolution integral equations of the form ##EQU1## where g(x,y) is the degraded image, f(x,y) is the desired ideal image, and p(x,y) is the known point spread function of the imaging system. The point spread function acts to blur or smooth out the ideal image, making it impossible to distinguish fine details in the recorded image g(x,y). Separately, g(x,y) is further contaminated by measurement noise. Thus: EQU g(x,y)=g.sub.e (x,y)+n(x,y), Equation (2)
where g.sub.e (x,y) is the blurred image that would have been recorded in the absence of noise, and n(x,y) represents the cumulative effects of all noise processes affecting final acquisition of the digitized array g(x,y). This includes the case of multiplicative noise, where n(x,y) is a function of f(x,y). The noise component n(x,y) is unknown but may be presumed small. Likewise, g.sub.e (x,y) is unknown. The type and intensity of the blurring caused by p, together with the magnitude of the noise in g, ultimately limit the quality of the restoration that can be achieved. It is convenient to write Equation (1) in operation notation as: EQU Pf=g Equation (3)
where P is the integral operation in L.sup.2 (R.sup.2) with kernel p(x-u,y-v).
The two dimensional Fourier transform plays a major role in the subsequent analysis. For a function a(x,y) of the space variables x,y, the Fourier transform a(.xi.,.eta.) may be expressed as: ##EQU2## The Fourier transform of the point spread function is called the optical transfer function and is an important tool in image analysis. The convolution theorem states that the transform of a convolution is the product of the individual transforms. The Parseval theorem states that the L.sup.2 scalar product of two functions in the space variables x,y is equal to their scalar product in the transform variables .xi.,.eta.. These and other results from Fourier analysis in L.sup.2 are well known in the art. For convolution equations such as Equation (1), it is advantageous to perform the analysis and computations in the Fourier domain. After image processing, the inverse Fourier transform is used to return to the space variables x,y.
The class G of point spread functions, described by Equation (5A) below, plays a key role in several civilian and military applications, including biomedical imaging; night vision systems; undersea imaging; imaging through the atmosphere; remote sensing; high definition television; and industrial applications. Consider first the class of point spread functions p(x,y) described as follows in the Fourier transform domain: EQU p(.xi.,.eta.)=e.sup.-.lambda.(.xi..spsp.2.sup.+.eta..spsp.2.sup.).beta., .lambda.&gt;0, 0&lt;.beta..ltoreq.1. Equation (5)
The case .beta.=1 corresponds to a Gaussian point spread function. This case occurs in quite diverse applications, including undersea imaging (see H. T. Yura, `Imaging in Clear Ocean Water,` Applied Optics, Vol. 12 (1973), pp. 1061-1066); low light-level electro-optical detection (see R. Weber, `The Ground-Based Electro-Optical Detection of Deep-Space Satellites,` Applications of Electronic Imaging Systems, Proceedings of the Society of Photo-Optical Instrumentation Engineers, Vol. 143, R. E. Franseen and D. K. Schroder, Eds. (1978), pp. 59-69); nuclear medicine gamma camera scintigrams (see S. Webbet al., `Constrained Deconvolution of SPECT Liver Tomograms by Direct Digital Image Restoration,` Medical Physics, Vol. 12 (1985), pp. 53-58; U. Raff et al., `Improvement of Lesion Detection in Scintigraphic Images by SVD Techniques for Resolution Recovery,` IEEE Transactions on Medical Imaging, Vol. MI-5 (1986), pp. 35-44; B. C. Penney et al., `Constrained Least Squares Restoration of Nuclear Medicine Images; Selecting the Coarseness Function,` Medical Physics, Vol. 14 (1987), pp. 849-859; B. C. Penney et al., `Relative Importance of the Error Sources in Wiener Restoration of Scintigrams,` IEEE Transactions on Medical Imaging, Vol. 9 (1990), pp. 60-70; K. S. Pentlow et al., `Quantitative Imaging of 1-124 Using Positron Emission Tomography with Applications to Radioimmunodiagnosis and Radioimmunotheray,` Medical Physics, Vol. 18 (1991), pp. 357-366); magnetic resonance imaging (see S. M. Mohapatra et al., `Transfer Function Measurement and Analysis for a Magnetic Resonance Imager,` Medical Physics, Vol. 18 (1991), pp. 1141-1144); and computed tomography scanners (see E. L. Nickoloff and R. Riley, `A Simplified Approach for Modulation Transfer Function Determinations in Computed Tomography,` Medical Physics, Vol. 12 (1985), pp. 437-442).
The case .beta.=5/6 describes blurring caused by atmospheric turbulence under long time exposure. This optical transfer function is important in post processing of degraded images obtained in airborne reconnaissance, in remote sensing from environmental or surveillance satellites, in the detection of deep-space satellites from ground-based telescopes, and in observational astronomy (see J. C. Wyant, Ed. Imaging Through the Atmosphere, Proceedings of the Society of Photo-Optical Instrumentation Engineers, Vol. 75 (1976); N. S. Kopeika, `Spectral Characteristics of Image Quality for Imaging Horizontally Through the Atmosphere,` Applied Optics, Vol. 16 (1977), pp. 2422-2426; R. Weber, `The Ground-Based Electro-Optical Detection of Deep-Space Satellites,` Applications of Electronic Imagining Systems, Proceedings of the Society of Photo-Optical Instrumentation Engineers, Vol. 143, R. E. Franseen and D. K. Schroder, Eds. (1978), pp. 59-69 N. S. Kopeika, `Imaging Through the Atmosphere for Airborne Reconnaissance,` Optical Engineering, Vol. 26 (1987), pp. 1146-1154; J. D. Gonglewski and D. G. Voelz, `Laboratory and Field Results in Low Light Post-detection Turbulence Compensation Using Self Referenced Speckle Holography,` Digital Image Synthesis and Inverse Optics, Proceedings of the Society of Photo-Optical Instrumentation Engineers, Vol. 1351, A. F. Gmitro, P. S. Idell, and I. J. LaHaie, Eds. (1990), pp. 798-806.
The case .beta.=1/2 corresponds to the Cauchy or Lorentzian density, and may be used to describe X-ray scattering in radiology (see F. C. Wagner et al., `A Characterization of the Scatter Point Spread Function in Terms of Air Gaps,` IEEE Transactions on Medical Imaging, Vol. 7 (1988), pp. 337-344).
A wide variety of electron-optical devices obey Equation (5) with a value of B satisfying 1/2.ltoreq..beta..ltoreq.1 (see C. B. Johnson, `A Method for Characterizing Electro-Optical Device Modulation Transfer Functions,` Photographic Science and Engineering, Vol. 14 (1970), pp. 413-415; C. B. Johnson, `Classification of Electron-Optical Device Modulation Transfer Function,` Advances in Electronics and Electron Physics, Vol. 33B (1972), pp. 579-584; C. B. Johnson et al., `High-Resolution Microchannel Plate Image Tube Development,` Electron Image Tubes and Image Intensifiers II, Proceedings of the Society of Photo-Optical Instrumentation Engineers, Vol. 1449, I. P. Csorba, Ed. (1991), pp. 2-12).
Such devices constitute important elements of various biomedical imaging modalities, including image intensifier-video camera (II-TV) fluoroscopic systems (see S. Rudin et al., `Improving Fluoroscopic Image Quality with Continuously Variable Zoom Magnification,` Medical Physics, Vol. 19 (1991), pp. 972-977); radiographic film digitizers (see F. F. Yin et al., `Measurement of the Presampling Transfer Function of Film Digitizers Using a Curve Fitting Technique,` Medical Physics, Vol. 17 (1990), pp. 962-966); radiographic selenium imaging plates (see P. J. Papin and H. K. Huang, `A Prototype Amorphous Selenium Imaging Plate System for Digital Radiography,` Medical Physics, Vol. 14 (1987), pp. 322-329); computed radiography photostimulable phosphor systems (see S. Sanada et al., `Comparison of Imaging Properties of a Computed Radiography System and Screen-Film Systems,` Medical Physics, Vol. 18 (1991), pp. 414-420; H. Fujita et al., `A Simple Method for Determining the Modulation Transfer Function in Digital Radiography,` IEEE Transactions on Medical Imaging, Vol. 11 (1992), pp. 34-39); digital TV tomography systems (see M. Takahashi et al., `Digital TV Tomography: Description and Physical Assessment,` Medical Physics, Vol. 17 (1990), pp. 681-685); and radiographic screen-film systems (see Y. Higashida et al., `Dual-Film Cassette Technique for Studying the Effect of Radiographic Image Quality on Diagnostic Accuracy,` Medical Physics, Vol. 11 (1984), pp. 646-652; H. Kuhn and W. Knupfer, `Imaging Characteristics of Different Mammographic Screens,` Medical Physics, Vol. 19 (1992), pp. 449-457).
Other important image tube/image intensifier applications include high definition television (HDTV) (see N. V. Rao, `Development of a High-Resolution Camera Tube for 2000-Line TV System,` Electron Tubes and Image Intensifiers, Proceedings of the Society for Photo-Optical Instrumentation Engineers, Vol. 1243, I. P. Csorba, Ed. (1990), pp. 81-86; R. Barden et al., `High Resolution MS-Type Saticon Pick-Up Tube with Optimized Electron-Optical Properties,` Electron Tubes and Image Intensifiers II, Proceedings of the Society for Photo-Optical Instrumentation Engineers, Vol. 1449, I. P. Csorba, Ed. (1991), pp. 136-147, and night vision and undersea imaging systems (see K. A. Costello et al., `Imaging GaAs Vacuum PhotoDiode with 40% Quantum Efficiency at 530 nm,` Electron Tubes and Image Intensifiers, Proceedings of the Society for Photo-Optical Instrumentation Engineers, Vol. 1243, I. P. Csorba, Ed. (1990), pp. 99-106.
In a typical imaging situation, several electron-optical devices may be cascaded together and used to image objects through a distorting medium such as the atmosphere or the ocean. The overall optical transfer function is then the product of finitely many functions of the type given by Equation (5), i.e., EQU p(.xi.,.eta.)=e.sup.-.SIGMA..spsp.J.sbsp.i=1.sup..lambda..sbsp.i.sup.(.xi.. spsp.2.sup.+.eta..spsp.2.sup.).spsp..beta.i, .lambda..sub.i.gtoreq. 0, 0&lt;.beta..sub.i .ltoreq.1. Equation (5A)
Problems having optical transfer functions satisfying Equation (5A) may be referred to as space-invariant image restoration problems wherein the image restoration problem is equivalent to an ill-posed backwards in time continuation problem for a generalized differential equation. Each factor in this product represents a component of the total blur, and is described by Equation (5) with a particular value for .lambda. and .beta.. A method for determining .lambda..sub.i and .beta..sub.i for each component is discussed in C. B. Johnson, `Classification of Electron-Optical Device Modulation Transfer Function,` Advances in Electronics and Electron Physics, Vol. 33B (1972), pp. 579-584. Equation (5A) also describes cascaded electron-optical devices forming the imaging chain in many digital radiography systems. Frequency domain characterization of such systems is an important and ongoing task, as is evident from the above-cited references. In several industrial applications, the general functional form defined by Equation (SA) may be found to best-fit an empirically determined optical transfer function, by suitable choices of the parameters .lambda..sub.i, .beta..sub.i and J.
It is emphasized that while many imaging phenomena do not have optical transfer functions that can be well-described by Equation (5A), the latter nevertheless encompasses a significant set of imaging problems. We denote by G the class of optical transfer functions defined by Equation (5A). Note that Equation (5) is a special case of Equation (5A). The functions in G are infinitely divisible probability density functions.
Most current approaches to image restoration are based on linear system inverse filter theory, where the input is the degraded image, and the single output is an approximation to the ideal image. Using a priori constraints on the unknown ideal image, various regularization methods are employed in these approaches, in order to stabilize the deconvolution procedure in the presence of noise. Tikhonov regularization is one of the best-known examples of such procedures. This method is closely related to Wiener filtering. The Tikhonov regularization method is also known as constrained least squares. In its simplest form, Tikhonov restoration is based on the following a priori information concerning the unknown ideal image f(x,y) in Equation (1), and the noise n(x,y) in the degraded image g(x,y): EQU .parallel.f.parallel..ltoreq.M, .parallel.n.parallel..ltoreq..epsilon., Equation (6)
where .parallel..parallel. denotes the L.sup.2 (R.sup.2) norm. This leads to the following constrained restoration problem for Equation (3). Given .epsilon., M&gt;0, with .epsilon./M&lt;&lt;1, find all f.epsilon.L.sup.2 (R.sup.2) such that EQU .parallel.Pf-g.parallel..ltoreq..epsilon., .parallel.f.parallel..ltoreq.M. Equation (7)
The a priori bound constraint on .parallel.f.parallel. in Equation (6) prevents explosive growth of noise in the restoration process. However, that constraint is not strong enough to prevent substantial noise contamination from obscuring fine detail in the restored image.
It is known in the prior art that the problem of restoring continuous dependence on the data in non well-posed problems for partial differential equations wherein the image restoration problem is equivalent to an ill-posed backwards in time problem for a generalized diffusion equation, requires constraints of the solution. This was taught by F. John, Continuous Dependence on Data for Solutions of Partial Differential Equations with a Prescribed Bound, Comm. Pure Appl. Math., 13 (1960), pp. 551-585. An effective approach in many cases is to prescribe an a priori bound M on the solution, or on the solution and finitely many of its derivatives, in an appropriate norm.
In the known prior art numerous authors have contributed to the study of stability estimates for improperly posed problems, linear and nonlinear, in various norms.
However, from the standpoint of numerical computation of the solutions from observed data, such estimates are not entirely satisfactory.
In the method of the present invention, a supplementary constraint is developed. This method is the slow evolution from the continuation boundary constraint, or sec B.
The theory of the SECB constraint is developed in detail in the manuscript Overcoming Holder Continuity in Ill-Posed Continuation Problems by Alfred S. Carasso, which is attached to this application and is an integral part of this application. The theory is discussed in the context of partial differential equations, and specifically in Sections 4, 4.1, 4.2, and 4.2 of the manuscript, with regard to image restoration. The basic algorithmic implementation of the SECB method, based on the Fast Fourier Transform, along with controlled computational experiments demonstrating the superiority of the SECB method over prior art, are discussed in Sections 4.2 and 4.3.
A research paper based on an expanded version of the above-mentioned manuscript was authored by Alfred S. Carasso and published in the SIAM Journal on Numerical Analysis, Volume 31, No. 6, December 1994, pp. 1535-1557.